require(cem)
require(survey)

load("analysis.RData")
#####MATCHING
analysis <- analysis[analysis$sample_resp_inp==1,]
analysis.match.1 <- subset(analysis, select=c(ohp_mo_treat, nnnnumhh_li_2,nnnnumhh_li_3,gender_inp,age_inp,dia_dx_pre_lottery_inp, hbp_dx_pre_lottery_inp,chl_dx_pre_lottery_inp,ami_dx_pre_lottery_inp,chf_dx_pre_lottery_inp,emp_dx_pre_lottery_inp,kid_dx_pre_lottery_inp,cancer_dx_pre_lottery_inp,dep_dx_pre_lottery_inp,edu_inp,hispanic_inp,race_white_inp,race_black_inp,race_nwother_inp,treatment,weight_total_inp,age_decile,bp_sar_inp,bp_dar_inp,bp_hyper,hbp_dx_post_lottery_inp,hbp_diure_med_inp, chl_inp, chl_high, hdl_inp, hdl_low, chl_dx_post_lottery_inp, antihyperlip_med_inp, a1c_inp, a1c_dia, dia_dx_post_lottery_inp, diabetes_med_inp, phqtot_high, dep_dx_post_lottery_inp, antidep_med_inp,cvd_risk_point_inp,any_dx_pre_lottery,age50_64,ohp_all_ever_inperson,rx_num_mod_inp,doc_num_incl_probe_inp,surg_num_incl_probe_inp,ed_num_incl_probe_inp,hosp_num_incl_probe_inp,chl_chk_inp,fobt_chk_inp,col_chk_inp,did_flu_inp,pap_chk_inp,mam50_chk_inp,psa_chk_inp,usual_clinic_inp,needmet_med_inp,med_qual_high_inp,smk_curr_bin_inp,obese,any_oop_spending_inp, tot_spend_other_trunc,catastrophic_exp_inp,owe_inp,borrow_inp,health_change_noworse,mcs8_score_inp,pcs8_score_inp,pain_low_inp,poshappiness_bin_inp, household_id))

### use these cutpoints and store as list 'cp' ####
cp<-list(nnnnumhh_li_2= c(0, 0.2, 0.4, 0.6, 0.8, 1),                                                                                                                                                                                                                                                                                      
nnnnumhh_li_3 = c(0, 0.5, 1),                                                                                                                                                                                                                                                                                                                 
gender_inp = c(0, 0.0588235294117647, 0.117647058823529, 0.176470588235294, 0.235294117647059, 0.294117647058824, 0.352941176470588, 0.411764705882353, 0.470588235294118, 0.529411764705882, 0.588235294117647, 0.647058823529412, 0.705882352941176, 0.764705882352941, 0.823529411764706, 0.882352941176471, 0.941176470588235, 1)  ,   
age_inp= c(20, 23.4, 26.8, 30.2, 33.6, 37, 40.4, 43.8, 47.2, 50.6, 54, 57.4, 60.8, 64.2, 67.6, 71)  ,                                                                                                                                                                                                                                  
dia_dx_pre_lottery_inp = c(0, 0.0714285714285714, 0.142857142857143, 0.214285714285714, 0.285714285714286, 0.357142857142857, 0.428571428571429, 0.5, 0.571428571428571, 0.642857142857143, 0.714285714285714, 0.785714285714286, 0.857142857142857, 0.928571428571428, 1)  ,                                                                          
hbp_dx_pre_lottery_inp = c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1) ,                                                                                                                                                                                                                                                                        
chl_dx_pre_lottery_inp = c(0, 0.0588235294117647, 0.117647058823529, 0.176470588235294, 0.235294117647059, 0.294117647058824, 0.352941176470588, 0.411764705882353, 0.470588235294118, 0.529411764705882, 0.588235294117647, 0.647058823529412, 0.705882352941176, 0.764705882352941, 0.823529411764706, 0.882352941176471, 0.941176470588235, 1) ,    
ami_dx_pre_lottery_inp = c(0, 0.0666666666666667, 0.133333333333333, 0.2, 0.266666666666667, 0.333333333333333, 0.4, 0.466666666666667, 0.533333333333333, 0.6, 0.666666666666667, 0.733333333333333, 0.8, 0.866666666666667, 0.933333333333333, 1) ,                                                                                                  
chf_dx_pre_lottery_inp= c(0, 0.0769230769230769, 0.153846153846154, 0.230769230769231, 0.307692307692308, 0.384615384615385, 0.461538461538462, 0.538461538461539, 0.615384615384615, 0.692307692307692, 0.769230769230769, 0.846153846153846, 0.923076923076923, 1)  ,                                                                               
emp_dx_pre_lottery_inp = c(0, 0.0588235294117647, 0.117647058823529, 0.176470588235294, 0.235294117647059, 0.294117647058824, 0.352941176470588, 0.411764705882353, 0.470588235294118, 0.529411764705882, 0.588235294117647, 0.647058823529412, 0.705882352941176, 0.764705882352941, 0.823529411764706, 0.882352941176471, 0.941176470588235, 1),     
kid_dx_pre_lottery_inp= c(0, 0.0769230769230769, 0.153846153846154, 0.230769230769231, 0.307692307692308, 0.384615384615385, 0.461538461538462, 0.538461538461539, 0.615384615384615, 0.692307692307692, 0.769230769230769, 0.846153846153846, 0.923076923076923, 1) ,                                                                                
cancer_dx_pre_lottery_inp= c(0, 0.0555555555555556, 0.111111111111111, 0.166666666666667, 0.222222222222222, 0.277777777777778, 0.333333333333333, 0.388888888888889, 0.444444444444444, 0.5, 0.555555555555556, 0.611111111111111, 0.666666666666667, 0.722222222222222, 0.777777777777778, 0.833333333333333, 0.888888888888889, 0.944444444444444, 1),
dep_dx_pre_lottery_inp=c(0, 0.142857142857143, 0.285714285714286, 0.428571428571429, 0.571428571428571, 0.714285714285714, 0.857142857142857, 1) ,                                                                                                                                                                                                   
edu_inp= c(1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4)   ,                                                                                                                                                                                                                                                          
hispanic_inp= c(0, 0.5, 1)   ,                                                                                                                                                                                                                                                                                                              
race_white_inp=c(0, 0.111111111111111, 0.222222222222222, 0.333333333333333, 0.444444444444444, 0.555555555555556, 0.666666666666667, 0.777777777777778, 0.888888888888889, 1) ,                                                                                                                                                             
race_black_inp=c(0, 0.25, 0.5, 0.75, 1) ,                                                                                                                                                                                                                                                                                                    
race_nwother_inp=c(0, 0.166666666666667, 0.333333333333333, 0.5, 0.666666666666667, 0.833333333333333, 1))


match<-cem(treatment="ohp_mo_treat",data=analysis.match.1, drop=c('household_id', 'ohp_all_ever_inperson','treatment', 'weight_total_inp', 'age_decile', 'bp_sar_inp', 'bp_dar_inp', 'bp_hyper', 'hbp_dx_post_lottery_inp', 'hbp_diure_med_inp', 'chl_inp', 'chl_high', 'hdl_inp', 'hdl_low', 'chl_dx_post_lottery_inp', 'antihyperlip_med_inp', 'a1c_inp', 'a1c_dia', 'dia_dx_post_lottery_inp', 'diabetes_med_inp', 'phqtot_high', 'dep_dx_post_lottery_inp', 'antidep_med_inp', 'cvd_risk_point_inp', 'any_dx_pre_lottery', 'age50_64', 'ohp_all_ever_inperson', 'rx_num_mod_inp', 'doc_num_incl_probe_inp', 'surg_num_incl_probe_inp', 'ed_num_incl_probe_inp', 'hosp_num_incl_probe_inp', 'chl_chk_inp', 'fobt_chk_inp', 'col_chk_inp', 'did_flu_inp', 'pap_chk_inp', 'mam50_chk_inp', 'psa_chk_inp', 'usual_clinic_inp', 'needmet_med_inp', 'med_qual_high_inp', 'smk_curr_bin_inp', 'obese', 'any_oop_spending_inp', 'tot_spend_other_trunc', 'catastrophic_exp_inp', 'owe_inp', 'borrow_inp', 'health_change_noworse', 'mcs8_score_inp', 'pcs8_score_inp', 'pain_low_inp', 'poshappiness_bin_inp'),eval.imbalance=T,cutpoints=cp)
match

write.csv(match$tab, file="match.table.csv")

md<-analysis.match.1[match$matched,]
md$weights<-with(md,match$w[match$matched])


#### Table 2 ####
svy.multi<-svydesign(~household_id, weights=md$weights, data=md)

T2.vars<-list("bp_sar_inp","bp_dar_inp", "bp_hyper", "hbp_dx_post_lottery_inp", "hbp_diure_med_inp", "chl_inp", "chl_high", "hdl_inp", "hdl_low", "chl_dx_post_lottery_inp", "antihyperlip_med_inp", "a1c_inp", "a1c_dia", "dia_dx_post_lottery_inp", "diabetes_med_inp", "phqtot_high", "dep_dx_post_lottery_inp", "antidep_med_inp","cvd_risk_point_inp")
T2.covs <- c("ohp_mo_treat+nnnnumhh_li_2+nnnnumhh_li_3+gender_inp+age_inp+dia_dx_pre_lottery_inp+ hbp_dx_pre_lottery_inp+chl_dx_pre_lottery_inp+ami_dx_pre_lottery_inp+chf_dx_pre_lottery_inp+emp_dx_pre_lottery_inp+kid_dx_pre_lottery_inp+cancer_dx_pre_lottery_inp+dep_dx_pre_lottery_inp+edu_inp+hispanic_inp+race_white_inp+race_black_inp+race_nwother_inp")
T2.form <- lapply(T2.vars, function(x) as.formula(paste0(x,"~", T2.covs)))
T2.models <- lapply(T2.form, function(x)  svyglm(x,design=svy.multi))
T2.eff <- lapply(T2.models, function(x) summary(x)$coef[2:6,1])
T2.se <- lapply(T2.models, function(x) summary(x)$coef[2:6,2])
T2.pval <- lapply(T2.models, function(x) summary(x)$coef[2:6,4])

T2 <-array(NA,dim=c(length(T2.vars),5,3),dimnames=list(c("BP Systolic", "BP Diastolic", "BP Elevated (%)", "Hypertension DX Post Lottery (%)","Current Hypertension Med Use (%)", "Total Cholesterol","High Cholesterol (%)", "HDL Level", "HDL Low(%)", "High Cholesterol DX Post Lottery (%)", "Current Use of Cholesterol Meds(%)","A1C Level","A1C Diabetic (%)","Diabetes DX Post Lottery (%)","Current Use of Diabetes Meds(%)","Depression Screening (%)","Depression DX Post Lottery","Current Use of Antidepressives (%)","CVD Risk Point"), c("1-6","7-12","13-18","19-24","25+"),c("Effect","SE","P-Value")))
T2[,,1]<-t(simplify2array(T2.eff))
T2[,,2]<-t(simplify2array(T2.se))
T2[,,3]<-t(simplify2array(T2.pval))
T2[c(3,4,5,7,9,10,11,13:16,18),,1]<-100*T2[c(3,4,5,7,9,10,11,13:16,18),,1]
T2[c(3,4,5,7,9,10,11,13:16,18),,2]<-100*T2[c(3,4,5,7,9,10,11,13:16,18),,2]

Table2<-matrix(NA,nrow=(dim(T2)[1]*dim(T2)[3]), ncol=(dim(T2)[2]+2))
r<-0
for (i in seq(1,3*dim(T2)[1],3)){
  r<-r+1
  Table2[i:(i+2),3:ncol(Table2)]<-rbind(T2[r,,1],T2[r,,2],T2[r,,3])
  Table2[i:(i+2),2]<-cbind(dimnames(T2)[[3]])
  Table2[i,1]<-dimnames(T2)[[1]][r]
}
colnames(Table2)[3:7]<-c(dimnames(T2)[[2]])


#write.csv(Table2,file="T2.multi.csv")


T3.vars<-list("health_change_noworse", "mcs8_score_inp", "pcs8_score_inp", "pain_low_inp", "poshappiness_bin_inp")
T3.covs <- c("ohp_mo_treat+nnnnumhh_li_2+nnnnumhh_li_3+gender_inp+age_inp+dia_dx_pre_lottery_inp+ hbp_dx_pre_lottery_inp+chl_dx_pre_lottery_inp+ami_dx_pre_lottery_inp+chf_dx_pre_lottery_inp+emp_dx_pre_lottery_inp+kid_dx_pre_lottery_inp+cancer_dx_pre_lottery_inp+dep_dx_pre_lottery_inp+edu_inp+hispanic_inp+race_white_inp+race_black_inp+race_nwother_inp")
T3.form <- lapply(T3.vars, function(x) as.formula(paste0(x,"~", T3.covs)))
T3.models <- lapply(T3.form, function(x) svyglm(x,design=svy.multi))
T3.eff <- lapply(T3.models, function(x) summary(x)$coef[2:6,1])
T3.se <- lapply(T3.models, function(x) summary(x)$coef[2:6,2])
T3.pval <- lapply(T3.models, function(x) summary(x)$coef[2:6,4])

T3 <-array(NA,dim=c(length(T3.vars),5,3),dimnames=list(c("Health same or better vs 1yr earlier (%)", "SF-8 Mental-component score","SF-8 Physical-component score", "No pain or very mild pain(%)", "Very happy or pretty happy(%)"), c("1-6","7-12","13-18","19-24","25+"),c("Effect","SE","P-Value")))
T3[,,1]<-t(simplify2array(T3.eff))
T3[,,2]<-t(simplify2array(T3.se))
T3[,,3]<-t(simplify2array(T3.pval))
Table3<-matrix(NA,nrow=(dim(T3)[1]*dim(T3)[3]), ncol=(dim(T3)[2]+2))
r<-0
for (i in seq(1,3*dim(T3)[1],3)){
  r<-r+1
  Table3[i:(i+2),3:ncol(Table3)]<-prettyNum(rbind(T3[r,,1],T3[r,,2],T3[r,,3]), digits=3)
  Table3[i:(i+2),2]<-cbind(dimnames(T3)[[3]])
  Table3[i,1]<-dimnames(T3)[[1]][r]
}
colnames(Table3)[3:7]<-c(dimnames(T3)[[2]])

#write.csv(Table3,file="T3.multi.csv")

T4.vars<-list("any_oop_spending_inp", "tot_spend_other_trunc","catastrophic_exp_inp", "owe_inp", "borrow_inp")
T4.covs <- c("ohp_all_ever_inperson+nnnnumhh_li_2+nnnnumhh_li_3+gender_inp+age_inp+dia_dx_pre_lottery_inp+ hbp_dx_pre_lottery_inp+chl_dx_pre_lottery_inp+ami_dx_pre_lottery_inp+chf_dx_pre_lottery_inp+emp_dx_pre_lottery_inp+kid_dx_pre_lottery_inp+cancer_dx_pre_lottery_inp+dep_dx_pre_lottery_inp+edu_inp+hispanic_inp+race_white_inp+race_black_inp+race_nwother_inp")
T4.form <- lapply(T4.vars, function(x) as.formula(paste0(x,"~", T4.covs)))
T4.models <- lapply(T4.form, function(x)  svyglm(x,design=svy.multi))
T4.eff <- lapply(T4.models, function(x) summary(x)$coef[2:6,1])
T4.se <- lapply(T4.models, function(x) summary(x)$coef[2:6,2])
T4.pval <- lapply(T4.models, function(x) summary(x)$coef[2:6,4])

T4 <-array(NA,dim=c(length(T4.vars),5,3),dimnames=list(c("Any out-of-pocket spending (%)", "Amount of out-of-pocket spending ($)","Catastrophic expenditures (%)", "Any medical debt (%)", "Borrowed money to pay bills or skipped payment (%)"), c("1-6","7-12","13-18","19-24","25+"),c("Effect","SE","P-Value")))

T4[,,1]<-t(simplify2array(T4.eff))
T4[,,2]<-t(simplify2array(T4.se))
T4[,,3]<-t(simplify2array(T4.pval))
Table4<-matrix(NA,nrow=(dim(T4)[1]*dim(T4)[3]), ncol=(dim(T4)[2]+2))
r<-0
for (i in seq(1,3*dim(T4)[1],3)){
  r<-r+1
  Table4[i:(i+2),3:ncol(Table4)]<-prettyNum(rbind(T4[r,,1],T4[r,,2],T4[r,,3]), digits=3)
  Table4[i:(i+2),2]<-cbind(dimnames(T4)[[3]])
  Table4[i,1]<-dimnames(T4)[[1]][r]
}
colnames(Table4)[3:7]<-c(dimnames(T4)[[2]])

#### Simulations ####

require(mvtnorm)
sim.models<-T2.models[c(1,2,3,4,5,6,7,10,11)]
names(sim.models)<-unlist(T2.vars)[c(1,2,3,4,5,6,7,10,11)]
var.cov<-lapply(sim.models,vcov)
sims<-vector("list",length(sim.models))
sims<-lapply(sim.models, function(x) rmvnorm(10000, na.omit(x$coef), sigma=vcov(x)))


sim<-lapply(sims, function(x) x[,2:6])
mean.sim<-lapply(sim,apply,2,mean)
quant.sim<-lapply(sim,apply,2,quantile, c(0.025,0.975))
in.sim<-lapply(sim,apply, 2, function(x) x[x>quantile(x,c(0.025)) & x<quantile(x,c(0.975))])

mean.sim<-lapply(mean.sim,function(x) c(0,x))
quant.sim<-lapply(quant.sim,function(x) cbind(0,x))
in.sim<-lapply(in.sim, function(x) cbind(0,x))
mean.sim<-lapply(mean.sim,function(x) 100*x)
quant.sim<-lapply(quant.sim, function(x) 100*x)
in.sim<-lapply(in.sim,function(x) 100*x)

png("simulate.figure1.png",width=800,height=400)
par(mfrow=c(1,2))
plot(c(1:6),mean.sim[[3]],type="o",ylim=c(-15,20), col="blue",pch=15, main="a. Hypertension Diagnosis", xlab="Months of Treatment", xaxt="n", ylab="Percent Difference")
legend("topleft",legend=c("% with hypertension", "% diagnosed with hypertension"),border=c("blue","green"),pch=c(15,16), col=c("blue","red"), bty="n")
polygon(c(rev(1:6),1:6),c(rev(quant.sim[[3]][1,]),quant.sim[[3]][2,]), border=NA, col="#0000ff22")
axis(side=1,at=1:6,labels=c("0","1-6","7-12","13-18","19-24","25+"))


lines(mean.sim[[4]], type="o", col="red", pch=16)
polygon(c(rev(1:6),1:6),c(rev(quant.sim[[4]][1,]),quant.sim[[4]][2,]),border=NA, col="#ff000022")

abline(h=0,lty=3,lwd=2)

plot(c(1:6),mean.sim[[3]],type="o",ylim=c(-15,20), col="blue",pch=15, main="b. Hypertension Medications",  xlab="Months of Treatment", xaxt="n", ylab="Percent Difference")
polygon(c(rev(1:6),1:6),c(rev(quant.sim[[3]][1,]),quant.sim[[3]][2,]), border=NA, col="#0000ff22")
legend("topleft",legend=c("% with hypertension", "% taking hypertension medication"),border=c("blue","dark green"),pch=c(15,17), col=c("blue","dark green"), bty="n")
axis(side=1,at=1:6,labels=c("0","1-6","7-12","13-18","19-24","25+"))

lines(mean.sim[[5]], type="o", col="dark green", pch=17)
polygon(c(rev(1:6),1:6),c(rev(quant.sim[[5]][1,]),quant.sim[[5]][2,]),border=NA, col="#11ff0029")

abline(h=0,lty=3,lwd=2)
par(mfrow=c(1,1))
dev.off()

png("simulate.figure2.png",width=800,height=400)
par(mfrow=c(1,2))
plot(c(1:6),mean.sim[[7]],type="o",ylim=c(-15,20), col="blue",pch=15, main="c. Cholesterol Diagnosis", xlab="Months of Treatment", xaxt="n", ylab="Percent Difference")
polygon(c(rev(1:6),1:6),c(rev(quant.sim[[7]][1,]),quant.sim[[7]][2,]), border=NA, col="#0000ff22")
axis(side=1,at=1:6,labels=c("0","1-6","7-12","13-18","19-24","25+"))
legend("topleft",legend=c("% with high cholesterol", "% diagnosed with high cholesterol"),border=c("blue","red"),pch=c(15,16), col=c("blue","red"), bty="n")

lines(mean.sim[[8]], type="o", col="red", pch=16)
polygon(c(rev(1:6),1:6),c(rev(quant.sim[[8]][1,]),quant.sim[[8]][2,]),border=NA, col="#ff000022")

abline(h=0,lty=3,lwd=2)

plot(c(1:6),mean.sim[[7]],type="o",ylim=c(-15,20), col="blue",pch=15, main="d. Cholesterol Medication", xlab="Months of Treatment", xaxt="n", ylab="Percent Difference")
polygon(c(rev(1:6),1:6),c(rev(quant.sim[[7]][1,]),quant.sim[[7]][2,]), border=NA, col="#0000ff22")
legend("topleft",legend=c("% with high cholesterol", "% taking cholesterol medication"),border=c("blue","dark green"),pch=c(15,17), col=c("blue","dark green"), bty="n")
axis(side=1,at=1:6,labels=c("0","1-6","7-12","13-18","19-24","25+"))

lines(mean.sim[[9]], type="o", col="dark green", pch=17)
polygon(c(rev(1:6),1:6),c(rev(quant.sim[[9]][1,]),quant.sim[[9]][2,]),border=NA, col="#11ff0029")

abline(h=0,lty=3,lwd=2)

par(mfrow=c(1,1))
dev.off()
